!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
program main
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_amr
   use interface_block_output
   use interface_check_state
   use interface_model_initial
   use interface_load_state
   use interface_save_state
   implicit none

   character(4) :: nt_char, cpu_char

   integer(4) :: ab_density, ab_pressure, i_level

   real(8) :: t_start, t_stop, t_wall
   real(8) :: max_speed_global

   type(block), pointer :: base, current

!----------------------------------------------------------------------------------------------------------------------------------|
!     mpi initialization
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_init(error)
   call mpi_comm_size(mpi_comm_world, cpu_number, error)
   call mpi_comm_rank(mpi_comm_world, cpu_id, error)
   cpu_id_p1 = cpu_id + 1

!----------------------------------------------------------------------------------------------------------------------------------|
   t_start = mpi_wtime()
   t_loop_e = t_start

!----------------------------------------------------------------------------------------------------------------------------------|
   if (cpu_id .eq. 0) then
      write (*, *) 
      write (*, '(a)') '|---------------|'
      write (*, '(a)') ' Program starts:'
      write (*, '(a)') '|---------------|'
      write (*, *)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   t = 0.0d0
   t_hc = 0.0d0
   t_wall = 0.0d0

!   open (unit_input, file = 'input')
!   read (unit_input, '(a)') directory
!   read (unit_input, '(a)') parameters_filename
!   close (unit_input)

   call getcwd(directory)
   directory = trim(directory) // '/'

!----------------------------------------------------------------------------------------------------------------------------------|
   t_s = mpi_wtime ()

   call read_parameters
   if (cpu_id .eq. 0) call system("mkdir -p " // trim(datamap_path))

   call mpi_barrier(mpi_comm_world, error)

!----------------------------------------------------------------------------------------------------------------------------------|
!   open (unit_datamap_manage_module, file = trim(directory) // '/data/data_manage_module.include')
!   write (unit_datamap_manage_module, '(a,i)') 'integer(4), parameter :: output_real8_flag = ', output_real8_flag
!   close (unit_datamap_manage_module)
   if (cpu_id .eq. 0) then
      open(unit_datamap_manage_module, file = 'data_manage_precision.include')
      write(unit_datamap_manage_module, '(a, i)') 'integer(4), parameter :: output_real8_flag = ', output_real8_flag
      close(unit_datamap_manage_module)
      open(unit_datamap_manage_module, file = 'data_manage_parameter.include')
      write(unit_datamap_manage_module, '(a, i, a, a, a)') 'character(', len(trim(parameters_filename)),                           &
         ') :: parameters_filename = ''', trim(parameters_filename), ''''
      write(unit_datamap_manage_module, '(a, i, a, a, a)') 'character(', len(trim(datamap_path)),                                  &
         ') :: datamap_path = ''', trim(datamap_path), ''''
      close(unit_datamap_manage_module)
   endif
!----------------------------------------------------------------------------------------------------------------------------------|
   allocate(base)
   allocate(base%x(nx_total), base%y(ny_total), base%child(ini_block_x, ini_block_y),                                              &
      base%p_pos(n_dims), base%start_pos(n_dims), base%end_pos(n_dims),                                                            &
      base%start_point(n_dims), base%id(n_dims, n_levels),                                                                         &
      base%fxro_bnd(ny_block, n_neighs), base%fxmx_bnd(ny_block, n_neighs),                                                        &
      base%fxmy_bnd(ny_block, n_neighs), base%fxmz_bnd(ny_block, n_neighs),                                                        &
      base%fxbx_bnd(ny_block, n_neighs), base%fxby_bnd(ny_block, n_neighs),                                                        &
      base%fxbz_bnd(ny_block, n_neighs), base%fxen_bnd(ny_block, n_neighs),                                                        &
      base%fyro_bnd(nx_block, n_neighs), base%fymx_bnd(nx_block, n_neighs),                                                        &
      base%fymy_bnd(nx_block, n_neighs), base%fymz_bnd(nx_block, n_neighs),                                                        &
      base%fybx_bnd(nx_block, n_neighs), base%fyby_bnd(nx_block, n_neighs),                                                        &
      base%fybz_bnd(nx_block, n_neighs), base%fyen_bnd(nx_block, n_neighs),                                                        &
      base%fxpo_bnd(ny_block, n_neighs), base%fypo_bnd(nx_block, n_neighs))

   call model_initial(base)

   call mpi_allreduce(max_speed, max_speed_global, 1, mpi_real8, mpi_max, mpi_comm_world, error)
   max_speed = max_speed_global

   base%dt = cfl / (max_speed / base%dx + max_speed / base%dy)
   base%dt = min(cfl * min(base%dx, base%dy) ** 2 / (max_et + eps), base%dt)
   base%dt_hc = min(cfl * min(base%dx, base%dy) ** 2 / (hc0 * max_te_ro + eps), base%dt) + eps
!   base%dt = min(cfl * min(base%dx, base%dy) ** 2 / (hc0 * max_te_ro + eps), base%dt) + eps
   min_dt = base%dt
   min_dt_hc = base%dt_hc

   if (load_state_flag .eq. 1) call load_state(base)

   t_e = mpi_wtime ()
   time_other = time_other + t_e - t_s

!----------------------------------------------------------------------------------------------------------------------------------|
   t_s = mpi_wtime ()

   if (cpu_id .eq. 0) then
      open (unit_par, file = trim(datamap_path) // '/params.mapout')
      write (unit_par, '(2x, a)') trim(directory)
      write (unit_par, '(i5, t64, a)') cpu_number, 'cpu_number'
      write (unit_par, '(i5, t64, a)') n_dims, 'n_dims'
      write (unit_par, '(i5, t64, a)') sys, 'sys'
      write (unit_par, '(i5, t64, a)') ng, 'ng'
      write (unit_par, '(i5, t64, a)') ini_block_x, 'ini_block_x'
      write (unit_par, '(i5, t64, a)') ini_block_y, 'ini_block_y'
      write (unit_par, '(i5, t64, a)') nx_no_bnd, 'nx_no_bnd'
      write (unit_par, '(i5, t64, a)') nx_block, 'nx'
      write (unit_par, '(i5, t64, a)') nx_total, 'nx_total'
      write (unit_par, '(i5, t64, a)') ny_no_bnd, 'ny_no_bnd'
      write (unit_par, '(i5, t64, a)') ny_block, 'ny'
      write (unit_par, '(i5, t64, a)') ny_total, 'ny_total'
      write (unit_par, '(i5, t64, a)') ratio, 'ratio'
      write (unit_par, '(i5, t64, a)') n_levels, 'n_levels'
      write (unit_par, '(i5, t64, a)') div_clean_flag, 'div_clean_flag'
      write (unit_par, '(f10.5, t64, a)') gm, 'gm'
      write (unit_par, '(f10.5, t64, a)') rangex, 'rangex'
      write (unit_par, '(f10.5, t64, a)') rangey, 'rangey'
      write (unit_par, '(f10.5, t64, a)') zero_posx, 'zero_posx'
      write (unit_par, '(f10.5, t64, a)') zero_posy, 'zero_posy'
      write (unit_par, '(f10.5, t64, a)') cfl, 'cfl'
      close (unit_par)
   endif

   t_e = mpi_wtime ()
   time_output = time_output + t_e - t_s

!----------------------------------------------------------------------------------------------------------------------------------|
   do while (.true.)

!----------------------------------------------------------------------------------------------------------------------------------|
!     remaining time
!----------------------------------------------------------------------------------------------------------------------------------|
      t_loop_e = mpi_wtime ()

      t_loop = t_loop_e - t_loop_s
!      if (ns .gt. 10 .and. cpu_id .eq. 0) then
!         if (t_loop .gt. t_loop_pre * 4.0d0) then
!            write (*, *)
!            write (*, *) 'Warning: expend too much time in one loop.'
!            write (*, *) 't_loop_pre: ', t_loop_pre
!            write (*, *) 't_loop: ', t_loop
!            write (*, *) 'Minimum number of blocks = ', minval (n_blocks)
!            write (*, *) 'Maximum number of blocks = ', maxval (n_blocks)
!            write (*, *) 'time_regrid one loop = ', time_regrid_one
!            write (*, *) 'time_integrate one loop = ', time_integrate_one
!            write (*, *) 'time_boundary one loop = ', time_boundary_one
!            write (*, *) 'time_output one loop = ', time_output_one
!            write (*, *) 'time_other one loop = ', time_other_one
!            write (*, *)
!         endif
!      endif
      t_loop_pre = t_loop
      t_remain = abs(t_end - t) / (min_dt / t_loop)
      t_loop_s = t_loop_e

      time_other = time_other + time_other_one
      time_output = time_output + time_output_one
      time_boundary = time_boundary + time_boundary_one
      time_regrid = time_regrid + time_regrid_one
      time_integrate = time_integrate + time_integrate_one
      time_regrid_one = 0.0d0
      time_integrate_one = 0.0d0
      time_boundary_one = 0.0d0
      time_output_one = 0.0d0
      time_other_one = 0.0d0

!      if (t .ge. nt * dt_output .or. ns .eq. 2) then
      if (t_output .ge. dt_output) then

!      call check_info ('output')

         t_s = mpi_wtime ()

         ab_density = 0
         ab_pressure = 0
         nt = nt + 1

         if (cpu_id .eq. 0) then
            write (*,"(a, i7, a, i6, a, es10.3, a, es10.3, a, es10.3)")                                                            &
               ' ns= ', ns, ', blocks= ', n_blocks_all, ', t= ', t, ', dt= ', min_dt
         endif

         write (nt_char, '(i4.4)') nt
         write (cpu_char, '(i4.4)') cpu_id

         if (cpu_id .eq. 0) then
            open (unit_nt, file = trim(datamap_path) // '/nt.mapout', form = 'unformatted')
            write (unit_nt) nt
            close (unit_nt)
            open (unit_t, file = trim(datamap_path) // '/info.t' // nt_char // '.mapout', form = 'unformatted')
            write (unit_t) t
            write (unit_t) min_dt
            close (unit_t)
         endif

         level_len = 0
         current => base%framework_next
         do while (current%lv .ne. 0)
            do i_level = 1, n_levels
               if (current%lv .eq. i_level) level_len(i_level) = level_len(i_level) + 1
            enddo
            current => current%framework_next
         enddo

         open(unit_datamap, file = trim(datamap_path) // '/data.t' // nt_char // '.cpu' // cpu_char // '.mapout',                  &
            form = 'unformatted')
         write(unit_datamap) n_blocks_vip_all
         write(unit_datamap) level_len
         write (unit_datamap) n_blocks(cpu_id_p1)

         current => base%next
         do while (current%lv .ne. 0)
            call block_output (current)
            current => current%next
         enddo
         close(unit_datamap)

         if (mod(nt, dn_save_state) .eq. 0 .and. cpu_id .eq. 0) then
            write (*, *)
            write (*, '(f6.2, a, i3, a, i3, a, f6.2, a)')                                                                          &
               t / t_end * 100d0, '% is completed, estimated remaining time:',                                                     &
               int(t_remain / 3600), 'h',                                                                                          &
               int(mod(t_remain, 3600.0d0) / 60), 'm',                                                                             &
               mod(t_remain, 60.0d0), 's'
            write (*, *)
         endif

         t_output = t_output - dt_output
         if (mod(nt, dn_save_state) .eq. 0 .and. save_state_flag .eq. 1) then
            n_save_state = n_save_state + 1
            call save_state(base)
            if (cpu_id .eq. 0) write (*,"(a, i6)") ' State save: ', n_save_state
         endif

         t_loop_pre = t_loop_pre * 100.0d0
         t_e = mpi_wtime ()
         time_output_one = time_output_one + t_e - t_s

         if (t .ge. t_end) exit

      endif

      max_et = 0.0d0
      call amr (base)

      t = t + min_dt
      t_output = t_output + min_dt

      t_s = mpi_wtime ()
!      call check_info ('other')
      max_speed = 0.0d0
      max_te_ro = 0.0d0
      current => base%next
      do while (current%lv .ne. 0)
         call check_state (current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en,         &
            current%x, current%y, nx_block, ny_block, gm, ab_density, ab_pressure, max_te_ro, max_speed)
         current => current%next
      enddo

      call mpi_allreduce (max_speed, max_speed_global, 1, mpi_real8, mpi_max, mpi_comm_world, error)
      max_speed = max_speed_global

!      base%dt = cfl / (max_speed / base%dx + max_speed / base%dy)
      base%dt = cfl * min(base%dx, base%dy) / max_speed / n_dims
      base%dt = min(cfl * min(base%dx, base%dy) ** 2 / (max_et + eps) / n_dims, base%dt)
      base%dt_hc = min(cfl * min(base%dx, base%dy) ** 2 / (hc0 * max_te_ro + eps) / n_dims, base%dt) + eps
!      base%dt = base%dt_hc
!      base%dt = min(cfl * min(base%dx, base%dy) ** 2 / (hc0 * max_te_ro + eps), base%dt) + eps

      base%t = base%t + min_dt

      t_e = mpi_wtime ()
      time_other_one = time_other_one + t_e - t_s

      ns = ns + 1
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  final step, free memory and output calculation time
!----------------------------------------------------------------------------------------------------------------------------------|
   if (associated (base)) deallocate (base)
   if (allocated (n_blocks)) deallocate (n_blocks)
   if (allocated (n_blocks_vip)) deallocate (n_blocks_vip)
   if (allocated (blocks_vip_flag)) deallocate (blocks_vip_flag)
   if (allocated (blocks_vip_flag_all)) deallocate (blocks_vip_flag_all)
   if (allocated (hilbert_start)) deallocate (hilbert_start)
   if (allocated (hilbert_end)) deallocate (hilbert_end)
   if (allocated (i_gene)) deallocate (i_gene)
   if (allocated (gene)) deallocate (gene)
   if (allocated (n_blocks_sendrecv)) deallocate (n_blocks_sendrecv)
   if (allocated (n_blocks_update_send)) deallocate (n_blocks_update_send)
   if (allocated (n_blocks_update_sendrecv)) deallocate (n_blocks_update_sendrecv)
   if (allocated (blocks_update_info_send)) deallocate (blocks_update_info_send)
   if (allocated (blocks_update_info_recv)) deallocate (blocks_update_info_recv)
   if (allocated (n_neighs_reflux_sendrecv_x)) deallocate (n_neighs_reflux_sendrecv_x)
   if (allocated (n_neighs_reflux_sendrecv_y)) deallocate (n_neighs_reflux_sendrecv_y)
   if (allocated (n_neighs_shareblock_sendrecv_x)) deallocate (n_neighs_shareblock_sendrecv_x)
   if (allocated (n_neighs_shareblock_sendrecv_y)) deallocate (n_neighs_shareblock_sendrecv_y)
   if (allocated (n_neighs_exchange_sendrecv_x)) deallocate (n_neighs_exchange_sendrecv_x)
   if (allocated (n_neighs_exchange_sendrecv_y)) deallocate (n_neighs_exchange_sendrecv_y)
   if (allocated (neighs_reflux_info_send_x)) deallocate (neighs_reflux_info_send_x)
   if (allocated (neighs_reflux_info_send_y)) deallocate (neighs_reflux_info_send_y)
   if (allocated (neighs_reflux_info_recv_x)) deallocate (neighs_reflux_info_recv_x)
   if (allocated (neighs_reflux_info_recv_y)) deallocate (neighs_reflux_info_recv_y)
   if (allocated (neighs_shareblock_info_send_x)) deallocate (neighs_shareblock_info_send_x)
   if (allocated (neighs_shareblock_info_send_y)) deallocate (neighs_shareblock_info_send_y)
   if (allocated (neighs_shareblock_info_recv_x)) deallocate (neighs_shareblock_info_recv_x)
   if (allocated (neighs_shareblock_info_recv_y)) deallocate (neighs_shareblock_info_recv_y)
   if (allocated (neighs_exchange_info_send_x)) deallocate (neighs_exchange_info_send_x)
   if (allocated (neighs_exchange_info_send_y)) deallocate (neighs_exchange_info_send_y)
   if (allocated (neighs_exchange_info_recv_x)) deallocate (neighs_exchange_info_recv_x)
   if (allocated (neighs_exchange_info_recv_y)) deallocate (neighs_exchange_info_recv_y)
   if (allocated (neighs_reflux_send_x)) deallocate (neighs_reflux_send_x)
   if (allocated (neighs_reflux_send_y)) deallocate (neighs_reflux_send_y)
   if (allocated (neighs_reflux_recv_x)) deallocate (neighs_reflux_recv_x)
   if (allocated (neighs_reflux_recv_y)) deallocate (neighs_reflux_recv_y)
   if (allocated (neighs_shareblock_send_x)) deallocate (neighs_shareblock_send_x)
   if (allocated (neighs_shareblock_send_y)) deallocate (neighs_shareblock_send_y)
   if (allocated (neighs_shareblock_recv_x)) deallocate (neighs_shareblock_recv_x)
   if (allocated (neighs_shareblock_recv_y)) deallocate (neighs_shareblock_recv_y)
   if (allocated (neighs_exchange_send_x)) deallocate (neighs_exchange_send_x)
   if (allocated (neighs_exchange_send_y)) deallocate (neighs_exchange_send_y)
   if (allocated (neighs_exchange_recv_x)) deallocate (neighs_exchange_recv_x)
   if (allocated (neighs_exchange_recv_y)) deallocate (neighs_exchange_recv_y)
   if (allocated (neighs_exchange_send_x_hc)) deallocate (neighs_shareblock_send_x_hc)
   if (allocated (neighs_exchange_send_y_hc)) deallocate (neighs_shareblock_send_y_hc)
   if (allocated (neighs_exchange_recv_x_hc)) deallocate (neighs_shareblock_recv_x_hc)
   if (allocated (neighs_exchange_recv_y_hc)) deallocate (neighs_shareblock_recv_y_hc)
   if (allocated (neighs_exchange_send_x_hc)) deallocate (neighs_exchange_send_x_hc)
   if (allocated (neighs_exchange_send_y_hc)) deallocate (neighs_exchange_send_y_hc)
   if (allocated (neighs_exchange_recv_x_hc)) deallocate (neighs_exchange_recv_x_hc)
   if (allocated (neighs_exchange_recv_y_hc)) deallocate (neighs_exchange_recv_y_hc)
   if (allocated (bnd_exterior)) deallocate (bnd_exterior)
   if (allocated (bnd_interior)) deallocate (bnd_interior)
   if (allocated (bnd_interior_pos)) deallocate (bnd_interior_pos)
   if (allocated (thr_va0)) deallocate (thr_va0)
   if (allocated (x_total)) deallocate (x_total)
   if (allocated (ro_total_x)) deallocate (ro_total_x)
   if (allocated (mx_total_x)) deallocate (mx_total_x)
   if (allocated (my_total_x)) deallocate (my_total_x)
   if (allocated (mz_total_x)) deallocate (mz_total_x)
   if (allocated (bx_total_x)) deallocate (bx_total_x)
   if (allocated (by_total_x)) deallocate (by_total_x)
   if (allocated (bz_total_x)) deallocate (bz_total_x)
   if (allocated (en_total_x)) deallocate (en_total_x)
   if (allocated (y_total)) deallocate (y_total)
   if (allocated (ro_total_y)) deallocate (ro_total_y)
   if (allocated (mx_total_y)) deallocate (mx_total_y)
   if (allocated (my_total_y)) deallocate (my_total_y)
   if (allocated (mz_total_y)) deallocate (mz_total_y)
   if (allocated (bx_total_y)) deallocate (bx_total_y)
   if (allocated (by_total_y)) deallocate (by_total_y)
   if (allocated (bz_total_y)) deallocate (bz_total_y)
   if (allocated (en_total_y)) deallocate (en_total_y)
   if (allocated (ro_total)) deallocate (ro_total)
   if (allocated (mx_total)) deallocate (mx_total)
   if (allocated (my_total)) deallocate (my_total)
   if (allocated (mz_total)) deallocate (mz_total)
   if (allocated (bx_total)) deallocate (bx_total)
   if (allocated (by_total)) deallocate (by_total)
   if (allocated (bz_total)) deallocate (bz_total)
   if (allocated (en_total)) deallocate (en_total)

   t_stop = mpi_wtime ()

   t_wall = t_stop - t_start

   if (cpu_id .eq. 0) then
      write (*, *)
      write (*, '(a)') '|-----------------|'
      write (*, '(a)') ' Program finished.'
      write (*, '(a)') '|-----------------|'
      write (*, *)
      write (*, '(a)') '|--------------------------------------------|'
      write (*, '(a, i3, a, i3, a, f6.2, a)') ' Total elapsed time:',                                                              &
         int(t_wall/3600.), 'h', int (mod (t_wall, 3600.) / 60), 'm', mod (t_wall, 60.), 's'
      write (*, '(a)') '|--------------------------------------------|'
      write (*, *)
      write (*, '(a)') '|--------------------------------------------------------------------------|'
      write (*, '(a)') ' Calculation costs:'
      write (*, '(a, f6.2, a)') ' Regrid:', time_regrid / t_wall * 100., '%'
      write (*, '(a, f6.2, a)') ' Integration:', time_integrate / t_wall * 100., '%'
      write (*, '(a, f6.2, a)') ' Data output: ', time_output / t_wall * 100., '%'
      write (*, '(a, f6.2, a)') ' Boundary condition: ', time_boundary / t_wall * 100., '%'
      write (*, '(a, f6.2, a)') ' Other (initial condition, data check, calculate minimum dt): ', time_other / t_wall * 100., '%'
      write (*, '(a)') '|--------------------------------------------------------------------------|'
      write (*, *)
   endif

   if (cpu_id .eq. 0) then
      open(unit_par, file = trim(datamap_path) // 'params.mapout', position = 'append')
      write(unit_par, '(a24, t32, i3, a, i3, a, f6.2, a)') ' Total elapsed time:',                                                      &
         int(t_wall/3600.), 'h', int (mod (t_wall, 3600.) / 60), 'm', mod (t_wall, 60.), 's'
      write(unit_par, '(a24, t32)') ' Calculation costs:'
      write(unit_par, '(a24, t32, f6.2, a)') ' Regrid:', time_regrid / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Integration:', time_integrate / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Data output:', time_output / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Boundary condition:', time_boundary / t_wall * 100., '%'
      write(unit_par, '(a24, t32, f6.2, a)') ' Other:', time_other / t_wall * 100., '%'
      close(unit_par)
   endif

   call mpi_finalize (error)

!----------------------------------------------------------------------------------------------------------------------------------|
   stop
end program main
